arXiv: 1509.02699v3 [astro-ph.SR] 4 Jan 2016 


Draft version January 5, 2016 

Preprint typeset using style emulateapj v. 5/2/11 


TOWARD MORE REALISTIC ANALYTIC MODELS OE THE HELIOTAIL: 
INCORPORATING MAGNETIC FLATTENING VIA DISTORTION FLOWS 

Jens Kleimann 

Ruhr-Universitat Bochum, Fakultat fiir Physik und Astronomic, Institut fiir Theoretische Physik IV, Bochum, Germany 


Christian Roken 

Universitat Regensburg, Fakultat fiir Mathematik, Regensburg, Germany 


Horst Fichtner 

Ruhr-Universitiit Bochum, Fakultat fiir Physik und Astronomic, Institut fiir Theoretische Physik IV, Bochum, Germany 


Jacob Heerikhuisen 

Department of Space Science and Center for Space Plasma and Aeronomic Research, 

University of Alabama in Huntsville, Huntsville, AL 35899, USA 
Draft version January 5, 2016 

ABSTRACT 

Both physical arguments and simulations of the global heliosphere indicate that the tailward 
heliopause is flattened considerably in the direction perpendicular to both the incoming flow and the 
large-scale interstellar magnetic field. Despite this fact, all of the existing global analytical models of 
the outer heliosheath’s magnetic field assume a circular cross section of the heliotail. To eliminate this 
inconsistency, we introduce a mathematical procedure by which any analytically or numerically given 
magnetic field can be deformed in such a way that the cross sections along the heliotail axis attain 
freely prescribed, spatially dependent values for their total area and aspect ratio. The distorting 
transformation of this method honors both the solenoidality condition and the stationary induction 
equation with respect to an accompanying flow field, provided that both constraints were already 
satisfied for the original magnetic and flow fields prior to the transformation. In order to obtain 
realistic values for the above parameters, we present the first quantitative analysis of the heliotail’s 
overall distortion as seen in state-of-the-art three-dimensional hybrid MHD-kinetic simulations. 

Subject headings: ISM: magnetic fields — magnetohydrodynamics (MHD) - methods: analytical - 
Sun: heliosphere 


1. INTRODUCTION 

Our knowledge about the large-scale structure of the 
heliosphere has recently been increased significantly with 
the Voyager 1 and 2 spacecraft that are by no w explor¬ 


ing in situ the local interstellar medium (LISM ; iGurnett 
et al.||2013l) and the inner heliosheath (see, e.g. Richard¬ 
son &: Burlaga 2013), respectively. This progress is ac- 


companied by ihsights gained from the remote measure¬ 
ments of energetic neutral atoms (ENAs) with the In- 
terstella r Boundary Exp l orer ( IBEX)', see the recent re¬ 
view by McComas et al. (2014). The ZBEAobservations 
supplernent those made with the Voyagers particularly 
regarding the global structure of the heliosphere because 
they are not limited to its upwind hemisphere but com¬ 
prise its entirety. While in this way the expected prin¬ 
cipal upwind-downwind asymmetry of the heliosphere 
could be confirmed, the actual structure of its down¬ 
wind hemisphere is still under debate. However, there 
can be no doubt that the dominant feature of the down¬ 
wind heliosphere is the so-called heliotail, along which 


jk@tp4.rub.de 

christian.roeken@mathematik.uni-regensburg.de 

h 1 @tp 4 .rub.de 

jacob.heerikhuisen@uah.edu 


the solar wind (SW) plasma eventually merges into the 
LISM. Conseq uently, t his he liotail — predicted to ex¬ 
ist already by Parker (1961) — is a natural feature of 
all classical models describing the interaction of the SW 
with the LISM. Nonetheless, its detailed structure has 
bee n the subje ct of only a few studies. In an early pa¬ 
per, Yu (1974) analyzed the magnetized “wake” of the 
SW under the influence of charge exchange with interstel¬ 
lar hydrogen atoms . This line of researc h wa s continued 
only much later by Jager & Fahr (1998) and Izmodenov 


& Alexashov (2003(1 While both papers conhrmed the 


significance of this process, the latter authors discussed 
a length of the heliotail of several tens of thousands AU, 
whereas the former fav ored a value of not more than 
1,500 AU. Additionally, I Jager & Fahr| (|1998|) recognized 
the potential significance of the heliotail tor the produc¬ 
tion of pick-up ions (PUIs). The rel a ted EN As were stud¬ 
ied by Czechowski & Grzedzielski| ( 19981), who further¬ 
more demonstrated how the direction and strength of an 
interstellar magnetic field (ISMF) can be deduced from 
the directi on of the heliota i l, an idea that was recently 
revived by McComas et al. (2013). Quantitative models 
of the forrn and deflection of the heliotail in the pres- 
ence of a magnetized LISM were first presented by |Ba-| 
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naszki ewicz fc Ratkiewicz ( 1989 ), M atsud a fc Fujimoto 
(|1993|) , and J Pogorelov &: Matsudal (I1998D. IKatkiewicz 
et ai.f (|2000|) , while not explicitly addressing the helio¬ 
tail topic, provided a classification of the asymmetries 
and distortions resulting from different ISMF orienta¬ 
tions and mentioned a flattening effect, which is also of 
relevance for the heliotail. Very recent studies looked at 
the stability of the heliotail (Po gorelov et al.|2014 2015|) 
and at the idea of its “splitting” ([Drake et al.|zU15 Opher 
et al. 2015), and, thereby, revived work on two top i cs tha t 
were also discussed to some extent already by|Yu|(|1974|). 

The heliotail is not only of significance for the tfux 
of PUIs and ENAs but also for that of cosmic rays 
(CRs). At the end of the last century, a so-called he- 
liomagnetotail anisotropy in the CR flux in t h e low - 
TeV range was discovered by Nagashima et al. (1998). 
This dipole-structured anisotropy and its — at least ap- 
proximate — relation to the heliotail was confirmed by 
the latest generation of large-area detectors like MILA- 
GRO, the Tibet Air Shower, IceCube, and others (Guil- 
lian et aLjpOO? Amenomori & Tibet Asy Gollaboration 
20101 |Karapetyan|2010[|Gngat fcT'ierre AUGER Gollab- 


oration 201 ip, showing an enhancement in the permille 

range of the low-TeV particle flux in the direction of the 
heliotail. A temporary denial of a physical link between 
the heliotail and this anisotropy was based on the facts 
(i) that the density in a gravitationally focused tail of in- 



magnetic reconnection as a process that accelerates CRs 
in the 50 GeV to 10 TeV range in the heliotail. While the 
existence of such a local source of CRs might be dou btful. 


in a subsequent paper, Desiati & Lazarian (2013) con¬ 
sidered the more likely anisotropy-inducing effect of the 
ISMF, whose homogeneity on the scale of th e heliosphere 


is dis turbed by the presence of the latter (see Roken et al. 


2015 and references therein). They claim that the large- 


scale CR anisotropy below 100 TeV is mostly shaped by 
particle interactions with turbulent ripples generated by 
the interaction of the heliospheric and the interstellar 
magnetic fields. 

None of these explanations or first “exploratory” mod¬ 
eling attempts were based on sophisticated heliospheric 
or advanced CR transport models, which are, however, 
both necessary in order to derive quantitative results 
that can be compared to observation s. This gap has 
recent ly be en fillede with the w ork by [Schwadron et ah 


(2014) and | Zhang et al. (2014), who studied the prob- 


lem m much more detail by parti cularly computing CR 


trajectories (as was first done by Washimi et al. 1999) 
rather than by employing the CR diffusion approxima- 
tion, which would represent a conceptual extreme in that 
case. These two papers arrived at the conclusion that a 
heliospheric impact on the CR anisotropy up to the TeV 
range must indeed be expected. Both studies, however, 
still have their drawbacks: Either an analytical model 


that is probably too simplifying ( Schwadron et al.l2014 ) 


or a numeri cal input, which is difhcult to handle ([Zhang 


et al. 2014), is used for the local ISMF configuration. 


For further analysis of the CR anisotropy problem, it is 
therefore desirable to have a significantly improved an¬ 
alytical model that incorporates crucial features of the 
ISMF. 

A first step to construct such a model was made by 


Roken et al. (2015), who derived an exact analytical so- 
lution for the local ISMF assumed to be frozen into the 
interstellar plasma flow. Not only this but, to the best 
of our kn owledge, all other analytical models of the lo- 
cal IS MF (|Whang|2dl0 Schwadron et al.|2014 Isenberg 


et al. 2015]) feature a circular cross section of the helio- 

tail, i.e. tar downtail, the heliopause takes the shape of a 
semi-i nhnite cylinder. None theless, the IBEX measure¬ 
ments (McC omas et al. 2013) as well as several numerical 
studies (e g. Heerikhuisen et al.||2014 Wood et al.||2014 | 
have confirrned the intuitive notion that with growing 
heliocentric distance, the heliotail becomes increasingly 
compressed perpendicular to the directions of the undis¬ 
turbed ISMF and the incoming interstellar flow. There¬ 
fore, we propose a method by which both the magnetic 
field B and a possibly associated flow held u, or in fact 
any other such vector held, may be deformed in a well- 
dehned manner into a realistic conhguration, such that 
it is still maintaining both the held line topology and 
the basic condition that the solenoidal magnetic held 
is frozen into the (possibly, but not necessarily equally 
solenoidal) plasma how. In other words, if 

( 1 ) 


and 


V • B = 0 


V X (u X B) = 0 


( 2 ) 


were satished by the model’s initial helds, then they will 
continue to be satished by the deformed vector helds. 

The outline of the paper is as follows. In Section 
the mathematical concept of distortion hows for the 
deformation of vector helds is introduced and the equa¬ 
tions relevant for the procedure are derived. Section 
provides some guidance for the choice of realistic param¬ 
eters by extracting the relevant data from self-consistent 
global simulations of the he liosphere, and applies them 
to the Roken et al. (2015) model. Finally, Section |4| 
completes the paper with a summary and conclusions. 
The fact that our method honors both the magnetic 
solenoidaliW condition ([^ and the stationary induction 
equation (|2| is proven rigorously in Appendix]^ while 
Appendices B] and [C] provide some details on the general 
distortion transformation and the employed data htting 
method, respectively. 


2. DERIVATION OF THE DISTORTION FLOW METHOD 

In order to deform a given pair of magnetic and 
velocity helds, we introduce the auxiliary held, wq, in 
which these helds are advected. The geometrical shapes 
of the distorted helds can be controlled by choosing a 
suitable form of wq. 


2.1. The Case of Constant Cross-seetional Areas 

In view of the intended application to the heliotail, we 
want to contract the initially circular cross-section in one 
direction and expand it in the perpendicular direction. 
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Figure 1. The flow fleld (l^, together with an initially circular 
set of points at time t = u (red) and the advected points at 
later time t = ti = 0.2/a > 0 (blue), when the distortion has 
reached an aspect ratio of 77 = exp(0.4) 1.49 according to Eq. 0- 

To this end, we fist consider the two-dimensional (2D) 
velocity field 


Wo : 




a X 

-ay 


( 3 ) 


with constant a > 0. The equation of motion for an 
inertialess particle being passively advected in this flow 

v{t) =wo[r(t)] , (4) 

the solution of which, subjected to the initial condition 
r(0) = ro, reads 


( \ _ ( xoexp{at) 

V ylt) ) ~ \ 2/oexp(-at) 


( 5 ) 


Specifically, an ensemble of particles which at t = 0 forms 
a unit circle will at t = > 0 have been deformed into 

an ellipse of half-axes 

(o,6) := (exp(ati),exp(—ati)) (6) 

with aspect ratio 

V ■=T = exp(2ati) > 1 (7) 

0 

while maintaining its area (oc a & = 1), in agreement with 
V • Wo = 0 (see Fig. [^. 

The idea is now to prescribe for each position z along 
the heliotail axis a “target aspect ratio” r]{z) = a{z)/b{z) 
into which the initially circular cross sections are to be 
deformed by an embedding of ([^ into three-dimensional 
(3D) space 

( 8 ) 
(9) 
( 10 ) 


as the solution to Eq. Q, evaluated at a fixed t = ti. 
This motivates a mapping —)■ 


X = a{z) xq 
y = b{z)yo 

Z = Zo 


( 11 ) 

( 12 ) 

(13) 


from arbitrary undeformed coordinates ro = {xo,yo,zo) 
at t = 0 to deformed coordinates r = (x, y, z) at t = ti, 
where a{z) and b{z) are given by Eq. (Iq with a = aiz). 
Note that this deformation has to be understood in the 
sense of a discrete transformation t = 0 —>■ ti, in co ntras t 
to the t-dependent, continuous case (cf. Section 2.3). 


Moreover, the requirement to ensure the constancy of 
the cross-sectional area tt a{z) b{z) for any z implies 


a(z) =77 (z) 1/2 and b{z) = 


(14) 


To see how a deformed vector P S {u,B} is obtained 
via the distortion field Wq, consider the unsqueezed vec¬ 
tor Pq located at rg and pointing in the direction of (5ro, 

i.e. Pg = vStq^ where Jrg is a small displacement vec¬ 
tor, which may also be pictured as a small segment of 
a field line. The constant v is introduced to warrant 
dimensional consistency, and can be used to accommo¬ 
date scalings of a desired magnitude. Via Eq. 0 , the X 
component of the deformed vector P at its new position 
r becomes 


Pxir) = iy[{x 5x) - a;] 

= ly [a(z Sz) (a;o -I- Sxq) — a(z) Xg] 
« ly [a(z) (jxg -I- a'(z) Xq 5 z\ 

= a{z) Poooiro) a'(z) Xq Poz(rg) , 


(15) 


where a prime denotes differentiation with respect to z, 
and the symbol has been used only because the 
terms of 0{S^) have been neglected in the Taylor expan¬ 
sion of a(z -I- Sz). The y component is found in complete 
analogy, and since the z component remains unaffected 
(i.e. zg = z), we end up with 


P. 


P. 


/ a Pi 
= bPr 


0 


Ox 

Oy 


Oz 


{a'/a) x 
{h'/b)y 


(16) 


In other words, given an arbitrary, unsqueezed vector 
field Pg, together with the function r]{z) specifying the 
desired aspect ratio relative to a circle, the components 
of the distorted field P at position r can be obtained by 


1. finding the corresponding 
rg(r) = (a;/o(z),y/5(z),z); 


“starting location” 


2. evaluating Pg(rg) at this position; and 


3. transforming the result according to Eq. (16). 


x{ti) = Xq exp [-l-a(zg) ti] 
y{ti) = yo exp [-a(zg) ti] 
z(ti) = zo 


For scalar quantities Q such as mass density or tem¬ 
perature, which contain no directional information, the 
proper form after distortion is simply Q{v) = (5o(rg). 
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Figure 2. Left: s elec t ed fi eld lines of the undistorted flux tube Middle: the corresponding magnetic field lines \2A\ after a distortion 

according to Eqs. | |19| l-| |20l l has been applied, viewed along the major (x) axis. Right: the same situation viewed from a different angle. 
The gray surface indicate s th e cylinder = 1 (left), which is transformed into the shape [x/a( 2 )]^ + \y/b{z)'^] = 1 via the 

squeezing transformation l|16[|. 


In the above derivation, the distortion flow wq (which 
mediates the deformation) is an auxiliary field that 
must be carefully distinguished from the actual velocity 
field u upon which the transformation acts. Wq is 
thus simply a mathematical object that acts similar 
to a flow field, although no physical flow of matter is 
associated with it. Likewise, the parameter t, which 
we have referred to as “time” for the sake of clarity, is 
rather to be understood as a deformation parameter 
controlling the extent of deformation from the original 
field Po S {uo,Bo} (at “time” t = 0) to its distorted 
counterpart P € {u,B} at t = ti. We continue to use 
the symbol t since the stationary induction equation ([^ 
implies that the configuration formed by u and B does 
not depend on physical time, such that no confusion 
should arise from this terminology. 


2.2. A Simple Example 

To illustrate the procedure, consider the cylindrical 
magnetic flux tube 


Bo(ro) 


Po 

1 + Po 




1 

1 + Po 


1 


-yo \ 

1 / 


(17) 


(where po := \/a;g + pg denotes the cylindrical radial co¬ 
ordinate and (po := arctan(?/g/a;o) the corresponding an¬ 
gular coordinate), together with a helical “swirl flow” 

/ - ft yo\ 

uo(rg) = flpo -b P , ( 18 ) 

V ^ / 

characterized by constants fl, P G K. It can easily be 
shown that these fields satisfy both Vq x (ug x Bg) =0 
and Vg • Bg = 0. In this example, we choose a deforma¬ 
tion with a spatially variable aspect ratio ri{z) = 1 + z^, 


which implies 


a(z) = \/i + z^ 
1 


b(z) = 


vT~ 


(19) 

( 20 ) 


according to Eq. (14). Since we are interested in the 
explicit forms of the distorted fields B and u, we now 
express all points (xo,yo,zo) in Pg(rg) through (x, i/, z) 
and apply the transformation (161 to the fields (17) and 
(flSl). Thus, with 


X y 

Px{x,y,z)=a{z) Pqx ( -yw, 


a(z) 


xPr 


a(z) ’ b(z) 

X y 


Oz 


a(z) ’ b(z )’ 


(21) 


, b'(z) 


biz] 


yPi 


a{z) ’ b{z )’ 


( X y 

Pz[x,y,z) = PQ^ —77-w,z 

Va(fl Kz) 

we obtain the distorted fields 
B= ^ 


\ + x"^ + y^(l -b z'^Y + 


u = 


1 + z^ 


( 22 ) 

(23) 


-y(l + Z^Y\ 


x-yz 

(24) 

1 + z^ ) 


/ xzV - y{l + zY"^ n\ 


xVl-yzV , 

(25) 

V {1 + zYV J 



which can ^ain be verified to satisfy the desired 
constraints ([^ and ([^. Fig. [2 illustrates the magnetic 
field configuration berore and mter the deformation. 
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2.3. Relation to Inductive Flux Transport 

The advection of field lines via distortion flows is rem¬ 
iniscent of the transport of frozen-in field lin es via the 
classical induct ion equation of ideal MHD (e.g.|Childress 
& Gilbert 1995). It is therefore reasonable to ask how the 
method detailed above relates to physical flux transport 
via an induction-type equation 


i9(P = V X (w X P) 


(26) 


that would lead to a similar result, however, for continu¬ 
ous deformations with respect to t described by the trans¬ 
formation (ro,t) !->■ (r{ro,t),t). To clarify this point, 
we first derive the explicit evolution equation that is ac¬ 
tually solved by the fields given by the distortion flow 
method. In the spirit of Eq. (15), but now with a vary¬ 
ing time/deformation parameter t, the vector field P(r, t) 
at position r = r(ro, t) and consecutive times t and t + 6t 
is 

P(r,t) =:/[r(r,t) - r(ro,t)] =vSL{r,t) (27) 
P(r,t -I- St) =v [r(r,t -I- St) - r(ro,t -|- (5t)] (28) 

= i/S'L{r,t + St) , 

where (5L denotes a small, time-dependent displacement 
vector, and f := r -|- dL. Note that t i—)■ r(ro,f) is the 
trajectory starting at rg, along which the P vector moves 
due to the action of the distortion flow w, yielding evolu¬ 
tionary deformation. Expanding (28) for small St, sub¬ 
tracting (27), and neglecting the terms of 0{St^), we 
obtain 


dtP = V\ WT - WT 


dr 

dt J~ ~ \dt dt ^ 

= V [w(r -I- (5L) — w(r)] 
« v^SP ■ V)w(r) , 
where the evolution equation 
dr 


(29) 


dt 


= w(r) 


(30) 


and a Taylor expansion of w(r -|- dL) for smal l up 
to first order have been used. Substituting Eq. (27) and 
again (30) into (29), one finds 


dtP = (P • V)w- (w- V)P 

= V X (w X P) — (V • P)w -)- (V • w)P , 


(31) 


which shows explicitly that the distortion flow method 
beco mes equivalent to solving the “induction equation” 
( 26 ) for the special case in which both V • w and V • P 


vanish. So in principle, one could as well use (26) to 
obtain some form of distorted fields, in which case 
9t(V • P) = 0 would even be satisfied unconditionally, 
i.e. without requiring any constraints on the value of 
V • w. However, as can be shown using an analysis 
similar to the one presented toward the concluding part 
of Appendix the alternative evolution equation (261 
can maintain V x (u x B) = 0 only for the special 
case V-P = 0 = V- w, in which it becomes identical 
to Eg. © anyway. Additionally, analytical solutions 
to (26) are usually difficult to obtain even for simple 


choices of w, whereas for the method proposed here, it 
suffices to solve the equation of motion (301, which will 


be much simpler to do in most cases. 


2.4. Generalization to Varying Cross-sectional Areas 

For a given application, the condition that cross sec¬ 
tions can change their shapes but not their absolute ar¬ 
eas may not always be satisfied. (This certainly holds 
true in the case of the heliospheric tail, as will be shown 
in Section 3.2 ) Therefore, we now wish to relax this 
constraint, such that a{z) and b{z) may be chosen inde¬ 
pendently of one another. Guided by our findings from 
the simpler situation considered so far, we will again take 
the distortion to be incompressible. This implies that the 
latter will no longer be confined to separate x-y planes, 
but that material from one z layer may be displaced into 
adjacent layers, i.e. the distortion flow corresponding to 
that of Eq. (§ will attain a non-zero z component. 

So consider a{z) and b{z) as given independently of 
one another. In an incompressible distortion, the spatial 
volume element is to be conserved, i.e. 


dxg dyo dzo = da; dy dz . 


(32) 


Additionally, we now allow for a prescribed displacement 
m{z) along the x axis as a third free parameter besides 
a{z) and b{z). The x and y components of the transfor¬ 
mation thus become 


Xo = 


— m(z) 
a(z) 


and 


2/0 = 


y 


b(z) ■ 


(33) 


A similar displacement in the y direction could easily 
be accounted for as well. We will refrain from doing so 
for simplicity of argument, and also because it will not 
be relevant for the intended application to the heliotail 
(see Section]^. While various choices of Zq = zo{x,y,z) 
would be consistent with Eq. (32), we restrict ourselves 
to the special case of zq = z^z), describing the sit¬ 
uation that all particles from the plane at zq are be¬ 
ing shifted into another plane at z. The general case 
of Zq = zo{x,y,z) is addressed in Appendix]^ Gonse- 
quently, the z component of the distortion transforma¬ 
tion has to satisfy the relation 


zo = 


d5 = 


i{S)b{S) di =: F{z) (34) 


in order to remain consistent with Eq. (32). As can be 


seen in this equation, z = 0 is to be interpreted as the 
“reference height” at which the transformation induces 
no vertical displacement of P, i.e. the planes zg = 0 and 
z = 0 coincide. 


Repeating the derivation of Eq. (15), we are led to 

Pa;(r) = IZ [(x d- Sx) — x] 

= V [a[z Sz) (xg -|- Sxq) d- m{z -\- Sz) 

- a{z) Xo - m{z)] (35) 

« iz [a{z) Sxo + a{z) Sz xq + m'{z) Sz] 

= a{z) Pox + [a{z) xg -I- m'{z)]P^ . 

The corresponding equation for Py is found in complete 
analogy, except that both m{z) and m'{z) are absent. 
Introducing now z = K(zo) as the inverse of zq = F{z), 
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we obtain 


P^(r) = iy[{z + dz) - z] 

= V [K{zo + 6zo) - K{zo)] 

= t^[dzoK{zo) Szo + 0(Sz'^)] 
« [a{z) b{z)]~^Po: 


(36) 


Oz 


for the z component. Inserting (36) into (35) and sup¬ 
pressing the argument z in a, b, and m, me complete 
vector transformation becomes 



Poz 


a Xq + m 

'’'r 


(37) 


Wh ile in this form it looks very similar to its analog 
(16), the most important difference is hidden in the 
additional requirement to obtain zq as a fun ctio n of z, 
which involves the evaluation of the in tegr al (34). This 
procedure will be illustrated in Section [372] 


3. REALISTIC PARAMETERS FOR THE HELIOTAIL 

Returning to our intended application, the heliotail, 
the above considerations naturally provoke the question 
of how the z profiles for the squeezing functions a, b 
and the displacement m should be chosen in order to 
obtain a configuration that resembles actual heliospheric 
conditions as closely as possible. 


3.1. Data Source and Fitting Method 

To provide some guidance for these choices, we turn 
to the 3D plasma-neutral simulations of the interac- 
tion between the S W an d the LISM performed by 
Heerikhuisen et al. ( 2014). These sim ulations employ 
an JVIHD solver (Pbgorelbv et al. 2006) for the ionized 
plasma and a Monte-Carlo part icle-based Boltzmann 
solver (Heerikhuisen et al. 2006) for neutral hydrogen. 
The neutral component is coupled to the plasma through 
charge-exchange collisions that generate energy and mo¬ 
mentum s ource terms which are applied to the MHD 
equations ( Heerikhuisen h Pogorelov|2010 ). For the sim¬ 
ulations considered here, the ion and neutral components 
were iterated until a steady-state solution to the SW- 
LISM boundary value problem was obtained. An impor¬ 
tant aspect of these simulations is the removal of energy 
from the SW through charge exchange with cold LISM 
hydrogen atoms as it flows down the heliotail. This leads 
to a reduction in pressure, and hence a decrease in the 
heliotail cross section at large distances. The simula¬ 
tions were performed on a spherical grid where the po¬ 
lar axis is aligned with the solar rotation axis, using an 
angular resolution of 3° in azimuth, 1.5° in colatitude, 
and a non-uniform radial grid with 284 cells. The inner 
boundary is located at 10 AU, and the outer boundary at 
1,000 AU. The resulting heliotail extends well beyond the 
outer boundary in the downstream LISM direction. The 
cell size increases with radial distance, which can make it 
difficult to accurately determine the heliopause location 
in the distant heliotail. However, the bulk properties 
of the plasma and neutral populations and the overall 
shape of the heliosphere are not significantly affected by 
this reduction in resolution near the outer boundary. 


Since these simulations were performed for magnetic 
field strengths of Hism G {1,2, 3,4} x 0.1 nT, they re¬ 
sulted in four distinct 3D configurations after a steady 
state had been reached. Here, these four datasets were 
then subjected to the following analysis: For all combi¬ 
nations j of magnetic field strengths i?ism and chosen x-y 
coordinate planes with z S {—100, —200,..., —800} AU, 
the heliopause locations Pj were determined as the iso¬ 
contours in the x-y plane at which the plasma tempera¬ 
ture has increased by a factor of five over its undisturbed 
LISM v alue at the outer comput ational boundary (see 
Fig. 6 in Heerikhuisen et al.||2014). Each Hj is then typ¬ 
ically given by iV « 200 pairs of coordinates {xi,yi)j, 
i = 1,..., iV, which have been rotated around the inflow 
(z) axis such that the x direction coincides with the di¬ 
rection of Bisin at infinity when projected onto the x-y 
plane. An ellipse 



(in which the second index j of Xi and yi has been sup¬ 
pressed), with major semi-axis Oj, minor semi-axis bj, 
and center coordinates (mj,m±j) was then fitted to each 
Hj contour by varying these four parameters until the 
total “difference area” Aj (i.e. the area covered by the 
interiors of either Hj or £j, but not both) attained its 
global minimum. The explicit form of Aj is given in 
Appendix [C} Here aj, bj, and ruj play the same role as 
the functions a, b, and m in the above derivation, except 
that they are no longer dimensionless but measured in 
absolute units of AU. Note that in order to obtain suf¬ 
ficiently good fitting results, we found it necessary also 
to allow for a displacement rn±j in the direction of f^’s 
minor axis as a fourth parameter. 

To illustrate this process, Fig. shows the fitting re¬ 
sults for the case of z = —500 AU) which was also used 


Heliopause shapes at z = —500 AU 



X [AU] 


Figure 3. Heliotail contours at z = —500 AU as derived from 
the hybrid simulation for -Bism values of 0.1 (red), 0.2 (blue), 
0.3 (green), and 0.4 nT (black), together with the respective 
fitting ellipses (dashed). The “+” symbol marks the inflow axis 
(x = 0 = y), while the “x” symbols indicate the centers of the 
respective ellipses. 
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Aspect ratio of heliotail [Heerikhuisen+ (2014) data] 



Figure 4. A spec t ratios rjj of fitting ellipses, together with 
fitting function ( |39[ l. 


Displacement along Bjg^^ [Heerikhuisen+ (2014) data] 



-z [AU] 


Figure 5. Linear displacement of the heliotail core axis against 
the inflow direction, together with fitting function l |40[ |. The 
dashed line marks a hypothetical deflection of 10° against the 
inflow axis. For the B = 0.4 nT case, data points at 2 < —4 00 AU 
were disregarded when setting up the fitting relation ( |40[ |. This 
is justified by the observation that for these cases the neTiopause 
contours develop a pronounced pear-shaped asymmetry that 
renders the elliptic approximation questionable. 


for Fig. 6 in Heerikhuisen et al. (20141. The minus sign 
honors the convention that the LiSM flow is inci dent 
from z —>■ + 00 , which will be exploited in Section |3.3[ 
The parameters rjj := aj/bj and rru thus derived are dis- 
played as functions of z in Figs. |^and[^ along with the 
heuristic fitting functions 


-2):= 1 + 3.27 


^looAuy 

1.5 


2.5 


z) := (487 AU) 


Si. 


nT 


(39) 


(40) 


which we propose to use as “best guesses” when setting 
up the deformation for a heliospheric field model. To our 
knowledge, this is the first time that simulations have 
been used to quantify the magnetic held-induced flatten¬ 
ing and deflection of the heliotail. 


While the cross-sectional flattening shows the expected 
behavior, i.e. it increases with both Bism and distance 
—z, the fact that the derived displacement does not de¬ 
pend monotonously on z is somewhat at odds with the 
notion of the heliotail axis being deflected by a fixed angle 
against the inflow axis (which would imply Wfit for 
fixed i?ism)- At this point, we may only note that while 
the simulations at hand do not support a correspond¬ 
ing scaling relation, the absolute magnitude of deflection 
is ind eed rather small, in agreement with [Wood et al.| 
(20141, who present observational constraints to argue 
for a deflection angle of the order of 10° and certainly 
below 20°. 

The resulting distribution of TO_L^fit(i?ism,-z) of abso¬ 
lute size (11 ± 6) AU and zero mean did not show any 
discernible systematics, and is therefore most likely due 
to the probabilistic nature of the hybrid MHD-kinetic 
approach (which can only accommodate a finite number 
of computational particles for the neutral component). 
Besides this, the only other cause for a departure of 
the problem’s symmetry with respect to the u-B plane 
at infinity would be the influence of the solar magnetic 
field, which, however, is apparently too weak to cause 
a systematic asymmetry of noticeable magnitude, in 
contrast to what has recently been claimed by |Opher| 


3 ' 

h; 


et al. (2015). Therefore, the perpendicular displacement 


to_l has not been considered further after fitting was 
completed. 


3.2. Cross-sectional Variations 

Since the distortion flow is to be incompressible, the 
constancy of the product a(z) 6(z) is a necessary criterion 
for the consistency of t he m ethod in its most basic form 
as presented in Section [O] It is therefore interesting to 
see to what degree the cross sections of our fitting ellipses 
can be considered constant. Fig. shows the effective 
cross-sectional radius Rj := ^/ojbjior the 0.3 nT case. 
The latter apparently increases toward a maximum at 
z « —500 AU and then decreases noticeably toward the 
far tail. This behavior may at least partly be attributed 


Effective radius of the heiiotaii, B|g,„ = 0.3 nT 



Figure 6. Variation in the heliotail’s effe ctiv e cross-sectional 
radius for the 0.3 nT case vs. fitting function 


























to the fact that the exact heliopause location is somewhat 
smeared in the simulations, especially in the distant tail, 
but also to variations of Uz and/or density along z, which 
are not captured in the simpler form of this approach of 
constant ab. Therefore, we proceed to ap ply t he more 
general deformation as outlined in Section |2.4[ thereby 
relaxing the requirement of cross-sectional constancy. 

To this end, we first note that according to Fig. the 
effective radius for the parameter range in question may 
be reasonably approximated by the heuristic function 


Rfit{z) ■■= Rc 



21 - 1/2 


(41) 


with best-fit parameters {Rc,Zc,Sc) = (314,-499,583) 
AU. This varying radius is obviously related to the di¬ 
mensionless squeezing functions a{z) and b(z) via 




(42) 


As can be seen from Fig. the cross sections can be ex¬ 
trapolated to an approximately circular shape at the ref¬ 
erence height 2: = 0, in agreement with a(0) = 6(0) = 1. 
This also means that it is indeed consistent to place the 
Sun in the z = 0 plane. From Eq. (34), we thus get 

[i?fit(0)]" zo= [ [Afit(5)]2 dz (43) 

Jo 


and hence 


zo 


1 + {Zc/bc)^ 

Zq — Cl 

with constants 


= Sr arctan 


arctan 


2: — 2:c 

2; — 2:c 


■Cn 


Co := arctan (zc/Sc) 

Cl := Sc + z^jSc 


(44) 

(45) 


(46) 

(47) 


which in this particular case evaluate to Co = 0.707 and 
Cl = 1,010 AU. If need be, this can immediately be 
inverted, yielding the explicit form oi z = K(zo) as 


3.3. Application to the Local ISMF 

As an illustrative application, we choose the rec ent an¬ 
alytical B fi eld model of the outer heliosheath by |Roken| 
et al. (2015), which provides the ex^t ISMF solution to 
the stationary induction equation (1^ for the case that 
the flow u derives from the Rankinenalf-body potential 
$ = uo{q/r + z) via u = —V$. Here, uq and dTruog, 
respectively, denote the LISM flow velocity at infinity 
and the SW source strength, and r := y/+ y'^ + z^ 
is the Sun-centered radial distance. The quantity ^/q 
may be interpreted as the upstream stand-off distance to 
the stagnation point, such that q = (125 AU)^ may be 
considered a reasonable choice for the solar case. The so¬ 
lution shows the expected behavior of the magnetic field 
piling up in front of and draping around the heliopause, 
which is easily identified as the surface satisfying 


H{r) := 2q — + y^) — Zy/Aq — (x^ + y'^) = 0 . (51) 

This model thus fulfills the criteria of ^ inclusion of both 
B and u fields which (ii) satisfy Eqs. ([^ and (|^, and (iii) 
feature a heliopause with circular cross section, and may 
therefore be subjected to the procedure outlined in the 
previous subsection. The distorted field B at position r 
is thus obtained by the following sequence of steps. 


1. Find the corres pon ding position rg in undistorted 
space via Eqs. (33) and (48), in which a(z), b(z ), 
and m{z) are fij^ using Eqs. (49)-(50) and (39)- 
(40). Since the upwind half-space z > 0 is taken 
tooe undistorted (implying ?7(z)|2>o = !)> the dis¬ 
torting transformation has to be limited to the tail 
region, i.e. the downwind half-s pace . Moreover, 
since the transformation formula (37) involves spa¬ 
tial gradients in a{z) and b{z) that would then lead 
to discontinuous B field components at z = 0, we 
introduce a spatial averaging function 


favg{z) ■— n 


1 -I- tanh 


(ra) 


(52) 


and replace all four functions /i(z) G {a(z), 6(z), 
to(z),E(z)} by 

y{z) := fi{z) [1 - favg{z)] + Ho /avg(z) , (53) 


z = Zc + Sc tan 



The two squeezing functions then become 


(48) 


a(z) 

6(z) 


^(06)^ = rfit(2;) y/ 77fit(0.3 nT, z) 
\ {a b)- = rot{z) 

V a -y?7fit(0.3 nT,z) 


(49) 

(50) 


such that all of the required information is available to 
turn any (heliospheric or other) magnetic field model 
with circularly constant cross sections into a more 
realistic version of itself by adjusting both the aspect 
ratio and the absolute cross-sectional area to a desired 
profile. 


thereby allowing them to tend smoothly to their 
undisturbed values (og, 60 , toq, Fb) = (1,1,0, zg), 
and ensuring that the distorted fields remain con¬ 
tinuous (and even differentiable) d ue to the now 
finite derivatives in transformation (37). 


2. The undistorted components of Bg are evaluated 
at this position in th is case using Eqs. (57)-(59) 
of Roken et al. (2015). 


3. The distorted components are computed from the 
undistorted ones using Eq. (37) in conjunction 
with the above parameters. Re mem ber that while 
to'(z) = 0 according to Eq. (40), the required 
derivatives a'(z) and b'{z) attain a more com¬ 
plicated form, which, however, may be obtained 
straightforwardly from their definitions. 
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Figure 7. Top row: perspective rendering of selected Bo(ro) field lines of the outer heliosheath model by |Roken et aL| ((2015|| for 
= 0.3 nT, as seen from three different vantage points. The gray surface is the heliopause as defined via H[ro) = U according to 
Eq. ( |51[ l, to which the magnetic field is tangential. The black lines on the heliopause are streamlines emanating from the stagnation point. 
Bottom row: the same situation for B(r) and H{r) = 0, i.e. after the distortion has been applied to both the magnetic field and the 
heliopause surface. Note that despite the massive change in shape, the magnetic field remains tangential to the heliopause. ||B|| values 
above 0.5 nT are capped. 


The resulting configuration can be seen in the side-by- 
side comparison of Fig. Evidently, the topological 
properties of the field have been maintained, while the 
shape of the heliopause has been deformed considerably. 


3.4. On the Possibility for Further Extensions 

In the deformations employed so far, the squeezing 
functions a and b have been allowed to vary in z but not 
within a given x—y plane. Therefore, by construction, 
they cannot accommodate asymmetric cross-sectional 
shapes which are found, e.g. in some of the simulation 
runs presented by Wood et al. (2014|. The method of 
distortion flows can, however, be easily generalized to 
such cases. To demonstrate this potential, we introduce 
a second parameter /3 > 0 and consider the distortion 
flow field E) to be replaced by the slightly more compli¬ 
cated field 


X 

y 




x^ 

-2xy 


(54) 


mapping —>■ K^, in which the y component of the /3 
term introduces a deformation that varies linearly in x 
(changing sign at the former symmetry plane a; = 0), 
and the x component of that term is chosen such that 
V • wi = 0 is satisfied. The corresponding, now coupled 
set of equations of motion r(t) = Wi[r(t)] has the ana¬ 
lytical solution 

x{t)=XQ\{xQl5+ l)e-x^{-at) - xqI5\ ^ (55) 

y(t) = yo exp(-Q;t) [xqP [exp(at) - 1] - l] ^ . (56) 

Evaluating these expressions at t = ti, and defining a 
new constant 

c:=(l/6-l)/3 (57) 

with b = exp(—ati) as before, we obtain the transforma¬ 
tion 

x = {xo/b){l - cxo)~^ (58) 

y={yob){l-cxo)'^ , (59) 

from which the geometric properties of the distorted con¬ 
tour may be easily derived. 


Wi : 
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Note, in particular, that while b may still be inter¬ 
preted as the contour’s absolute extension along the y 
axis in both directions, the position of maximum exten¬ 
sion is shifted away from {x, y) = ( 0 , ± 6 ) toward positive 
X, while in the perpendicular direction, the contour ex¬ 
tends to ([(c±l) 6]“^,0) rather than (±&“^,0). This 
also implies that a spatially bounded and simply con¬ 
nected contour requires c < 1. Hence, b continues to 
adjust the scaling oc ( 1 / 6 , 6 ), while the parameter c con¬ 
trols the degree of deformation, ranging from the usual 
ellipse (c = 0 ) to something like a rounded triangle for 
values near unity. 

After some algebra, the resulting vector transformation 
turns out to be 


^ (l + bcx)^ ^ , {b^c'x-b')x ^ 

- 7 ^Qx \ 7 J 02 

b b 

Py = 2c{l + bcx)y Po. + Poy 

-I- [l -I- ( 6 c -I- 6 'c -I- 26c')a;] y Pqz 


(60) 

(61) 


P. = Po. (62) 

since z = zq when embedding the 2D flow field ( [^ 
into 3D space. The proof of explicit compliance with 
constraints Q and may be carried out in close 
analogy to the one presented in Appendix]^ and is not 
given here. Fig. illustrates how the parameters 6 and 
c may be used to approximate an asymmetric heliopause. 



-500 0 

X (AU) 



500 -500 0 500 

X(AU) 



Figure 8. Upper panel: cross-sectional shapes of the heliopause 
a.t z = —500 A U (left) and z = —2 000 AU (right), as seen in the 
simulations by |Wood et al.| ||2014|l, adapted from Fig. 5 of that 
work. Lower panel: two distorted contours (blue) as approxi¬ 
mations to the simulated cross sections, plotted relative to the 
reference circle corresponding to 6 = 1 and c = 0 (red), which has 
the same total area due to V • wi = 0. The respective parameters 
employed for the left and right contour are (6 = 0.85, c = 0.10) 
and {b = 0.60, c = 0.15). 


4. SUMMARY AND CONCLUSIONS 

We have introduced and described the mathematical 
procedure of distortion flows, by which any pair of 3D 
magnetic and velocity fields B and u can be deformed 
into different geometric shapes — but with the same 
topological properties — that are more desirable in a 
given application. It was proven rigorously that the 
method conserves both V x (u x B) = 0 and V • B = 0. 

For our principal application, namely, the heliospheric 
tail, we employed sophisticated global simulations to de¬ 
rive heuristic profiles for the former’s cross-sectional as¬ 
pect ratio and deflection as functions of position along 
the tail axis, and used this information to deform a recent 
analytic tail model into a still analytic, yet more realistic 
version of itself using the distortion flow method. 

It should, however, be kept in mind that vector equa¬ 
tions different from Eqs. Q and (|^ are not guaranteed to 
remain unaff ected by the squeezing t ransformation. For 
instance, the Schwadron et al.| (|2014|) heliosphere model 
is built around the requirement that there be no currents 
outside the heliopause (V x B = 0), and while our ap¬ 
proach could certainly be used to deform their magnetic 
field, the resulting B field will most likely have VxB 7 ^ 0. 

We note that the use of a (solenoidal) distortion flow 
does not imply any constraints on the solenoidality (or 
lack thereof) of the physical fields that are being dis¬ 
torted. In particular, we do not in any way argue for the 
plasma to be incompressible. 

Furthermore, we would like to stress that the appli¬ 
cability of the distortion flow approach outlined in this 
paper is suitable for any source magnetic field, and does 
not depend on its particular physical context. Therefore, 
it is a rather generic tool that could clearly be applied 
outside the heliospheric context as well. In view of the 
rather small number of known analytic solutions to the 
MHD induction equation, the method clearly has the 
potential to yield the solutions for the magnetized flow 
around complex geometries simply by starting from a 
known solution of s impler geometry (like a sphere, see, 
e.g. the appendix of Isenberg et al. (2015) and references 
therein) and then deforming this solution as desired. 


ACKNOWLEDGMENTS 

We are grateful to Gunnar Hornig, Yuri Litvi¬ 
nenko, Frederic Effenberger, and Andreas Kopp for 
helpful discussions. We acknowledge financial support 
via the project El 706/15-1 funded by the Deutsche 
Forschungsgemeinschaft (DFG). J.H. acknowledges sup¬ 
port from NASA grants NNX12AB30G, NNX14AF43G, 
NNX14AF43G, NSF grant OGI-1144120, and Depart¬ 
ment of Energy grant SC0008334. We also appreciate 
discussions at the workshop “Reconnection, Turbulence, 
and Particles in the Heliosphere” in Queenstown, New 
Zealand, as well as at the team meeting “Heliosheath 
Processes and Structure of the Heliopause: Modeling 
Energetic Particles, Cosmic Rays, and Magnetic Eields” 
supported by the International Space Science Institute 
(ISSI) in Bern, Switzerland. 


































11 


APPENDIX 


A. INVARIANCE OF SOLENOIDALITY AND FROZEN-IN CONDITIONS 


We proceed to prove that if the vector fields u and B are computed using transformation (37) from fields Ug and 
Bq satisfying Vq • Bq = 0 and Vq x (uq x Bq) = 0, these two equations in deformed coordinates will also hold for u 
and B. To this end, we first need to see how the differential operator V = {dx,dy,dz) in distorted space is related to 
its undistorted counterpart Vg = {dxo,dyg,dzg). From Eqs. (33) and ([M|, we have 


1/a 0 0 \ 

0 1/6 0 

— {a'x + m')/a —{yb')/b ab J 


-'Vo 


(A.1) 


which is found by explicit calculation via 


dx dx dxo dx dyo dx Ozq dx \ a{z) ) dxo a(z) dxo 

d dxo d dyo d dzo d n ^ ^ \ ^ n 

dy dy dxo dy dyo ^ dy dzo ^ dy V6(z) j dyo ^ b{z) dyo 
d dxo d dyo d dzo d 

dz dz dxo dz dyo dz dzo 

d /x —m(z)\ d d f y \ d dF{z) d 

dz \ a{z) J dxo dz V6(z) / dyo dz dzo 

{x — m{z)) a'{z) + m'{z) a{z) d yb'{z) d / \ \ 9 

=- m ^ 

_ xa'{z) + m'{z) d yb'{z) d ^ ^ ^ ^ 

a{z) dxo b{z) dyo dzo - 


Then, by means of transformation (37), we obtain 


„ dPx dPy 

V • P = ' ^ 


dz 


a Pf 


Ox 


dx dy 
l_d_ 
a dxo 

a' xo + m' d 
a 


a' Xo + m' 
a b 

h'yo d 


Pqz 
+ a 6 


dPpx J_ 
dxo a b 
1 

a b 


dxo b dyo ' dzp 

d fa' xo + m' \ d 
dxo \ a + dyo 

a' Xo + m' dPoz \ , f b'yo dPpz 


1 d 
b dyo 
d \ Pqz 
a b 


bPoy + 


b'yo 
a b 


naz 


P 


= Vg • Pg + ^ 


dxo J \ b 

d f a' Xo + m' 

dxo V ® 

'-V- 

—a'ja 


dyo 
d 


+ 


dyo 


b'yo 


Pf 


Oz 


dP 


Oy 


dyo 


+ 


dP 


Oz 


dzo 


■ Poz a b . 

dzo \ab 


b'yo 


-P 


Oz 


1 da 1_^ 
a dzo b dzo 


=b'/b 


in which the final pair of square brackets cancels the preceding one due to 


a 




da{z) 

dz 


a b 


da 

dzo 


(A.2) 

(A.3) 


(A.4) 


(A.5) 


(A.6) 


and similarly for b' = dzb{z). This proves the invariance of the solenoidality condition for P. 

For the induction equation, we start by defining the undistorted electric field Eg := —Ug x Bg, such that Vg x Eg = 0 
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according to Eq. ([^, and then evaluate 

—E = u X B 


auox + {a'xo + m')/{ah) \ /a Bq^ + {a'xo + m')/{ah) \ 

buoy + {b'yo)/{ab) mqz x b Boy + {b'yo)/{ab) Bq^ 

uo^/{ab) J \ Bo:,/{ab) J 

(uoyBoz - uozBoy)/a 
{"^OzBox "^OxBqz)/ b 

(uoxBoy - uoyBox) ab - (itoy^Oz - uozBoy)ia'xo + m')/a - {uozBqx - uqxBqz) b'y^/b 
( Bqx IOj \ 

Eoy/b 

\ Eo^ab- Eox (a'xo + rn')la - Eoy b'yo/b / 


and further 
V X E = I 


(l/a) dxo 

i^/b) dy„ 


\ -(a'xo + rn')/a dxo “ (b'yo)/b dyo+ab dzg 
The X component reads 


(V X E), = 


1 _d_ 
b dyo 


a'xo + m' ^ 
txQz — i^OfC- 


Oy 


b'yo 


/ Eox/a 

X Eoy/b 

V Eoz ab- Eox (a'xo + m')/a - Eoy b'yo/b . 


a'xo + m' d b'yo d , d \ Eoy 


dEoz dEoy \ f dEoy dEox ^ a'xo + m' f a db b' 


dyo dzo J \ dxo dyo 


a b 


b dzo 62 


~ — I Eoy — 0 


— (VoXEo)cc 


—(VoxEo)2 


(A.7) 


(A.8) 


(A.9) 


in which the last term vanishes for the same reason as above. The y component is computed in complete analogy, and 
vanishes as well. Finally, 


1 d f Eoy\ 1 d /Ei 


a dxo 




1 / dEoy dEox 


b dyo \ a J ab\ dxo dyo 


= 0 


(A.IO) 


— (VoxEo)z 


which concludes the proof of the invariance of the induction equation under transformation (37). 


The above considerations can be extended to more general distortion flow fields w as follows. Since in this case the 
transition from Pq to P proceeds via the evolution equation (31), it is convenient for the proof of the invariance of 
constraint ([^ to consider 


dt(S/ • P) = V • (atP) = V • [V X (w X P) - (V • P)w + (V • w)P] 

= 0 - [(V • P)(V • w) + w V(V • P)] + [(V • 'W)(V • P) + P • V(V • w)] (A.ll) 

= P-V(V-w)-w-V(V-P) . 


If V • w = 0, then this is the well-known advection partial differential equation (PDE) 

{dt+w-V)f = 0 (A.12) 

for f(r,t) := V • P, meaning that / is passively advected (and thus constant) along flow lines of w. In particular, 
f(x,t) = 0 is the only solution that satisfies the appropriate initial condition /(r, 0) = (V • P)|t=o = 0 for all r. This 
shows explicitly that V • P continues to vanish as long as V • w = 0. In fact, we see that any w with constant (not 
necessarily vanishing) divergence, like for instance w cx r, will conserve the solenoidality of P. This is intuitively 
plausible, since such a linear scaling corresponds to an isotropic “magnification” of space that leaves all angles between 
vectors, etc. unchanged. 

Regarding the second constraint 0, we first define the vector quantity to be conserved as C := V x (u x B), for 
which we need to show that for the initial condition C|i=o = Oj 


dtC = at(V X (u X B)) = V X ([9tu] xB + ux [9tB]) 

= V X ([(u • V)’w — ("w • V)u] X B -f u X [(B • V)w — (w • V)B]) 
vanishes at any t > 0. Using the identities 
[(A.V)B]^= A,dyB, 

X (A X B)j ^ = Cijk dj(SklmAlBxn/} — (bilSjm bimbjl) dj(AiBzn) — dj(AiBj AjBj/) 


(A.13) 


(A.14) 

(A.15) 
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for general vectors A, B S (with the convention £123 = 1, and summation over double indices is implied), we obtain 
for the ith component of Eq. ( |A.13[ ) 

OtCi ^k^k^j^^i\ 

“t“ ^Uii^Bj^dkWj WkdkBj'^ Uj(^Bj^dkWi 'Wk^kBi^^^ 

— “b i^iBk 'fkkB^^i^dkkUj^ dki^UjBj^ U^Bj^Wk^ 

= djiukBj - UjBk){dkWi) + {ukBj - UjBk){djkWi) +dj{utBk - UkB^){dkWj) 


(A.16) 


[Vx(uxB)]fc 

“b {kiiBk UkBi^ (^djkWj^ 
[V(V-w)];, 


= :S -.N 

- djkiujBi - UiBj) Wk + dk{ujBi - UiBj){djWk) 


-afe[Vx(uxB)]i 


-N 


in which the second ter m S va nishes due to the symmetry of djk and the summation over both indices. Upon reverting 
to vector notation, Eq. (A.16) thus simply becomes 


dtC = (C • V)w- (w • V)C 


(A.17) 


when evaluated for the special case V • w = 0 that was already fo und ne cessary for the invariance of the solenoidality 
constraint. Since w does not depend on t, one can easily use Eq. (A.17) to show via mathematical induction that 


(arC)|t=o = 0 VneNo. 


(A.18) 


The PDE system (A.17), together with the initial conditions (A.18), states a Cauchy problem. Under the real ana- 
lyticity assumption on all coefhcients of this PDE system, or rather the distortion vector field w, one can apply the 
Cauchy-Kovalevskaya theorem for local existence and uniqueness. From this, it follows that there exists a unique local 
analytic solution of the Cauchy pr oblem , namely, for C, in the neighborhood of t = 0. Consequently, the Taylor series 
of C at t = 0 with the conditions (A.18) indeed yields the trivial solution 


cw = E 


n—O 


1 a"c 

^ dt^ 


n\ 


t" = 0 


(A.19) 


t=o 


proving that if C is zero initially, it stays zero at all later times. Again, a vanishing gradient of V • w is sufficient 
to obtain this result. This is analogous to the usual MHD induction equation for 9tB, which evolves a pre-existing 
magnetic field B non-linearly but cannot do so without an initial seed field. 

We not e in pa ssing that it would not be permissible to conclude already from C|(=o = 0 and the fact that the linear 
equation (A.17) admits the trivial solution C = 0 that C has to vanish iden tically, as is sometim es done erroneously 
in the literature for the usual, time-dependent MHD induction equation (e.g. Subramanian||2000 ). This can clearly be 
seen by the simple example of the ordinary linear differential equation g{t) = ( 2 /t) g[t)^ which possesses the general so¬ 
lution g{t) = gi that is non-zero for t > 0 and any constant gi = ^(l) G IR\{0}, despite the fact that g{0) = 0 = 5 ( 0 ). 


B. THE CASE OF ARBITRARY VERTICAL DISPLACEMENTS 


The PDE for the determination of the new 0 component of the general coordinate transformation defined by Eq. (33) 
and Zq = Zo(x, y, z), satisfying volume conservation, can be obtained by first computing the exterior derivatives of me 
original coordinates xq, yo, and Zq, reading 


dsn = - 


dyo = 


dx — 


[x — to] a' 


-b to' I dz 


dy - 


yb' 


d^o = ^dx 
ox 


b 

dzp 

dy 


dz 

dy 


dzo 

dz 


dz 


The corresponding volume three-form becomes 


dxo A dyo A dzo = ^ 
a b 


[a; — to] a' 


dzp 

dx 


yb' dzp 
b dy 


dzp 

dz 


da; A dy A dz 


in which the symbol A denotes the standard alternating wedge product on the exterio r alg ebra 
volume conservation, i.e. da;o A dyo A dzp = da; A dy A dz, we immediately find from Eq. (B.2) that 


[x — m]a' 


dx b ay 5z “ “ ■ 


(B.l) 


(B.2) 

Requiring 

(B.3) 
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This PDE is the defining equ ation for the general component transformation zn = z n(x, y, z). In the special case 
zo = zo(z) discussed in Section 2.4 we have dzo/dx = 0 = dzojdy, and thus from (B.3) we obtain 


dzp 

dz 


= ab . 


(B.4) 


Integration leads to 


zo = / a(z) b{z) dz , 
Jo 


(B.5) 


which is the component transformation used in Eq. (34). 


Let 


C. ELLIPSE FITTING PROCEDURE 


(fiij := arctan 


/ - mj_j \ 
\ Xi- ruj ) 


(C.l) 


be the angle between the x axis and the vector which points from the center of Ej to the contour point (xi, S 'Kj. 
The respective distances from the center to that point and to the intersection with in the same direction are then 
given by 


= [{xi - rrijf + {yi - m±jf] 

^ X 2 / . ^ 21 “1/2 

COS (fij 




/ sin 


(C.2) 

(C.3) 


Using the identities cos^ = 1/(1 + tan^ y}) and sin^ = (tan^ (/?)/(! + tan^ yE) allows us to write the quantity to be 
minimized as 


N 

i-1 


N 


2=1 


1 



(C.4) 


The contour points {xi,yi)j are generated with equidistant distribution of arctan(j/i/a;i)j. Strictly speaking, A_, would 
need 5(pij := <Pi+ij — to be uniform for all i in order to actually be proportional to the total difference area 
between T-Lj and £j. However, given that there is in any case some freedom in the details of a chosen fitting method, 
this discrepancy is certainly small enough to be safely ignored for the present purpose. 
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